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ABSTRACT 



The trailing vortices generated by the control planes of sub- 
marines give rise to surface signatures in the form of scars and 
striations. 

Two counter-rotating vortices were generated in a novel experi- 
mental system and their interaction with the free surface was investi- 
gated. In addition, the governing equations have been solved through 
the use of the boundary-element method for a representative Froude 
number. The results have been expressed in terms of the depth of 
submergence of the vortices, their mutual induction velocity, and the 
initial vortex separation. It has been shown that the free surface 
begins to deform when the vortices are at a distance of about one ini- 
tial vortex separation from the free surface. The height of the maxi- 
mum deformation is attained at a normalized time of about 0.1, when 
the vortices are at a distance of about 0.5 b D from the free surface. 
The elevated part of the surface is bounded by two scars, whose 
motion is slaved to that of the vortices. 



TABLE OF CONTENTS 



I. INTRODUCTION 13 

II. PHYSICAL ANALYSIS 17 

A DIMENSIONAL ANALYSIS 17 

B GOVERNING DIFFERENTIAL EQUATIONS AND 

BOUNDARY CONDITIONS 19 

C. DIMENSIONLESS PARAMETERS 24 

1. Dynamic Scale 24 

2. Buoyant Scale 28 

III. NUMERICAL SOLUTION TECHNIQUES 34 

A INTRODUCTION 34 

B FINITE DIFFERENCE TECHNIQUES 35 

C "FINITE ELEMENT METHODS 41 

D. BOUNDARY INTEGRAL EQUATION METHODS 44 

E. HYBRID METHODS 47 

F. NUMERICAL METHODS. CONCLUSIONS, AND 

SELECTION , 49 

IV. DISTRIBUTED VORTEX MODEL 54 

A INTRODUCTION 54 

B. APPLICATION OF THE DISTRIBUTED VORTEX METHOD . 56 



5 



C RESULTS OF THE NUMERICAL CALCULATIONS 



59 



V. EXPERIMENTS 61 

A EXPERIMENTAL APPARATUS 61 

R PROCEDURES 61 

VI. DISCUSSION OF RESULTS 63 

VII. CONCLUSIONS 65 

APPENDIX: FIGURES 66 

LIST OF REFERENCES. 108 

INITIAL DISTRIBUTION LIST 112 



6 



LIST OF FIGURES 



1 Computational Domain, Governing Equations, 
and Boundary Conditions 



2.1 


Velocity 


Field 


for 


T* 


= 


-2.40 


2.2 


Velocity 


Field 


for 


T* 


= 


-2.30 


2.3 


Velocity 


Field 


for 


T * 


= 


-2.20 


2.4 


Velocity 


Field 


for 


T* 


= 


-2.10 


2.5 


Velocity 


Field 


for 


T * 


= 


-2.00 


2.6 


Velocity 


Field 


for 


T* 


= 


-1.90 


2.7 


Velocity 


Field 


for 


T* 


= 


-1.80 


2.8 


Velocity 


Field 


for 


T* 


= 


-1.70 


2.9 


Velocity 


Field 


for 


T* 


= 


-1.60 


2.10 


Velocity 


Field 


for 


T* 


- 


-1.50 


2.11 


Velocity 


Field 


for 


T* 


= 


-1.40 


2.12 


Velocity 


Field 


for 


T* 


= 


-1.30 


2.13 


Velocity 


Field 


for 


T* 


= 


-1.20 


2.14 


Velocity 


Field 


for 


T* 


= 


-1.10 


2.15 


Velocity 


Field 


for 


T* 


= 


-1.00 


2.16 


Velocity 


Field 


for 


T* 


= 


-0.90 


2.17 


Velocity 


Field 


for 


T* 


= 


-0.80 


2.18 


Velocity 


Field 


for 


T* 


= 


-0.70 


2.19 


Velocity 


Field 


for 


T* 


= 


-0.60 


2.20 


Velocity 


Field 


for 


T* 


= 


-0.50 


2.21 


Velocity 


Field 


for 


T* 


= 


-0.40 



66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 
87 



7 



2.22 Velocity Field for T* = -0.30 88 

2.23 Velocity Field for T* = -0.20 89 

2.24 Velocity Field for T* = -0.10 90 

2.25 Velocity Field for T* = 0.1 0 91 

2.26 Velocity Field for T* = 0.20 92 

2.27 Velocity Field for T* = 0.10 93 

3 Experimental Apparatus 94 

4. 1 Vortex Motion at T* = -4.05 95 

4.2 Vortex Motion atT* = -3.46 96 

4.3 Vortex Motion at T* = -3.00 97 

4.4 Vortex Motion atT* = -2.73 98 

4 . 5 Vortex Motion at T* = -1 .63 99 

4 . 6 Vortex Motion at T* = -0.82 100 

4.7 Vortex Motion at T* = -0.28 101 

4 . 8 Vortex Motion at T* = -0.09 102 

4.9 Vortex Motion at T* = 0.09 103 

4.10 Vortex Motion at T* = 0.28 104 

4.11 Vortex Motion atT* = 0.44 105 

4.12 Vortex Motion at T* = 0.82 106 

5 Free Surface Shape at Maximum Rise 107 



8 



TABLE OF SYMBOLS AND ABBREVIATIONS 



b 0 Initial vortex pair spacing 

d Q Initial depth of the vortex pair 

F(t) .Arbitrary function of time 

F v Froude number 

g Gravitational acceleration 

i vn 

Lc Characteristic length 

Ln Integral scale of the turbulent field 

n N2/N? 

N c Characteristic Brunt-Vaisala frequency 

N 0 Brunt-Vaisala frequency 

P Pressure 

P a -Ambient Pressure 

q Scalar velocity Vu 2 + v 2 

qm q/uc 

r e Effective core radius of the vortex 

Re Reynolds number 

r m Nondimensional radial distance 

Ri, R 2 Radii of curvature of the free surface 
t Time 

t m Nondimensional time (u c t/L c or tN c ) 

u X-component of velocity 



9 



U c Characteristic velocity 

u m u/U c 

v Y-component of velocity 

v m v/U c 

V 0 Initial mutual induction velocity 

w(z) .Complex velocity function 

x Horizontal component of the coordinate axis (parallel to the 

line joining the vortex cores) 

x m x/Lc 



y Vertical component of the coordinate axis 

ym y/Lc 

x’ m , y' m Position of vorticity 
zi Complex position of a vortex 

z 0 Complex position of the right-hand vortex in the vortex pair 

Zi Complex conjugate of zi 

z 0 Complex conjugate of z 0 

[B] Forcing function matrix 

[K] Coefficient matrix 

[U] Matrix of unknowns 

e Rate of decay of turbulent energy per unit mass 

e* Nondimensional e 

T m Circulation of a given distributed vortex 
T o Circulation of the vortex pair 

rj Position of the free surface 

Tim n/Lc 



10 



p Density of water 

Pm p/ Po 



p 0 Reference density of water in stratified medium 

p(y) Initial variation of density with y 

p'(x,y,t) Fluctuating part of density with time 
0 -Velocity potential 

0 m 0/r o 



K 

V 

c 

Cm 

V 2 



Geometric pie 

Kinematic viscosity of water 

Vorticity 

CLc/Uc 

fd 2 d 2 \ 

Laplacian Operator I^T + ^rl 

Nondimensional Laplacian operator 





a Surface tension 



1 1 



ACKNOWLEDGMENT 



The author wishes to express his sincere thanks to Distinguished 
Professor T. Sarpkaya for his invaluable assistance, guidance, and 
patience throughout this investigation and especially during the many 
months of program development. I’m sure to him it seems more like 
years. It has been a true privilege to share some of this man’s genius 
and understanding of nature. I consider this opportunity to be the 
hallmark of my experience at the Naval Postgraduate School. 

Also, the author wishes to thank Mr. Jack McKay of the Mechan- 
ical Engineering Department Machine Shop for his hard work and 
unique craftsmanship in building and operating the experimental 
apparatus. 



12 



I. INTRODUCTION 



Vortices and vortex wakes have become a major theme of 
aerodynamics research since the advent of the large aircraft and the 
understanding of their evolution required an examination of many of 
the fundamental problems in fluid mechanics. Much of the progress 
made during the past two decades was discussed at the Symposium on 
Aircraft Wake Turbulence and Its Detection [Ref. 1] and the Aircraft 
Wake Vortices Conference [Ref. 2]. Comprehensive reviews of the 
entire subject have been given by Donaldson and Bilanin [Ref. 3], 
Widnall [Ref. 4], and Hallock and Eberle [Ref. 5]. 

These studies, as well as numerous others carried out since 1977, 
have uncovered a number of complex problems which must be 
resolved in order to achieve a better understanding of the important 
features of trailing vortices in homogeneous and stratified media. The 
principal ones are as follows [Ref. 6]: 

1. Roll-up process: The velocity and turbulence distribution at any 
station behind the wing depend on the wing section, wing-tip 
shape, Reynolds number, wing incidence, and the distance of 
the station from the wing [Ref. 7]. The distributions of the initial 
velocity and turbulence, which influence the roll-up and the 
decay process, cannot be changed independently. For example, 
a change in the tip shape changes the core size, as well as the 
velocity and turbulence distributions. High levels of turbulence 
result in an increased diffusion of vorticity, which in turn 
increase the core size. 

2. Probe sensitivity of the vortices: Flow visualization studies 

suggest that trailing vortices are extremely sensitive to 
disturbances created by even very small probes or bubbles. This 
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forces one to use non-intrusive means of measurement such as a 
Laser Doppler Velocimeter. Even then, “vortex wandering” 
(Ref. 8], which makes the vortices appear larger than normal in 
time-averaged velocity measurements (for vortices generated by 
a wing in a wind tunnel), or the unsteady nature of the flow (for 
vortices generated by a wing in a tow basin) makes the mean 
velocity profiles in the vortices difficult to determine. 

3. Large-scale instabilities: The vortices are seldom observed to 
decay away owing to viscous and turbulent dissipation, but are 
almost always destroyed by either mutual induction instability 
(Crow instability (Ref. 9]) and/or vortex breakdown. The Crow 
instability grows exponentially, and results either in linking of 
the vortex pair into a series of crude vortex rings or in a highly 
disorganized intermingling of the vortices. 

Vortex breakdown, whose mathematical details have not yet 
been adequately treated, rearranges the vortex structure and 
increases the core size, turbulence, and energy dissipation. 
Thus, it is very difficult to measure accurately the trajectories of 
the three-dimensional vortices from their creation to their 
ultimate demise. 

4. Reynolds number: Even the highest Reynolds numbers, based . 
on wing chord, reached in wind tunnels or towing basins, are an 
order of magnitude lower than what is possible for an aircraft. 
Thus, the scale effects are not easy to assess. 

5. Ambient conditions such as turbulence and stratification play 
major roles in the evolution of vortices. The quantification of 
these effects requires numerical analysis and extremely careful 
experiments (Ref. 10]. 

6. Ground or free surface effects: The vortex pair may move toward 
a rigid boundary at which the no-slip condition must be satisfied 
or toward a free surface at which the zero-shear condition must 
be satisfied. In either case, the vortices come under the 
influence of their images and move accordingly. 

The phenomenon is further complicated by several additional 
facts. When the vortices are propelled toward a rigid surface, vorticity 
of opposite sign is generated on the no-slip boundary and swept 
toward the vortex pair. The total vorticity diminishes very quickly as 
vorticity from the two regions diffuses, the wall region serving as a 
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strong sink for the vorticity associated with the original vortex 
[Ref. 11]. The development of a boundary layer along the rigid wall 
may give rise to flow separation for sufficiently high Reynolds numbers. 
With or without such a separation, howfever, the center of the vortex 
pair eventually moves away, or “rebounds,” from the wall [Refs. 11- 
13]. 

For the case of a zero-stress boundary, the free surface still acts as 
a vorticity sink, but this is relatively weak due to the absence of 
intense oppositely signed vorticity. Thus, in the absence of other 
impending phenomena, one expects a mild interaction between the 
vortices and the free surface and a small rebound of the vortex pair 
from the free surface. However, the ability of the free surface to 
deform under the influence of strain fields leads to a strong inter- 
action between the vortices and the free surface. 

It is evident from the foregoing that the motion and the life-span 
of trailing vortices are governed by a number of nonlinearly dependent 
complex phenomena. A number a experimental and analytical studies 
have been carried out at the Naval Postgraduate School by Sarpkaya 
and his students [Refs. 6, 14-21] in order to investigate the effects of 
these parameters on the rise and demise of the trailing vortices in 
homogeneous and density stratified media. These studies have clearly 
identified the various demise mechanisms in both media and estab- 
lished basic relationships between the rise of vortices and the 
governing parameters in a finite as well as effectively infinite medium 
[Ref. 20], free from ambient turbulence. 
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The present investigation is a continuation and refinement of the 
previous studies. The intent is to analyze the rise and demise of a 
vortex pair in a medium with a deformable free surface. This problem 
is of interest both from the standpoint of determining the interaction 
effect of the free surface on the rise and demise of the vortex pair and 
from the standpoint of predicting the resulting free surface shape- 
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II. PHYSICAL ANALYSIS 



A. DIMENSIONAL ANALYSIS 

The dependent parameter of major importance to the problem 
solution 'is the instantaneous position of the vortex pair (x, y). It may 
be expressed as a function of the following parameters [Ref. 22]: 

x = f(t, V 0 , d 0 . p 0 . dp/dy, v, b Q , g. r e , e. Ln) (1) 



and 

y = f(t, V 0 , d 0 . p 0 . dp/dy, v, b G , g. r e . e, Ln) (2) 

in which the variable definitions are as follows: 
t time 

V 0 initial mutual induction velocity of the vortices 

d 0 “initial depth of the vortex pair 

Po reference density of the medium 

dp/dy linear density gradient 

v kinematic viscosity of the medium 

b 0 initial separation of the vortex pair 

g gravitational acceleration 

r e effective core radius of the vortex 

e rate of decay of the turbulent energy per unit mass 

Ln ' integral scale of the turbulent field 

The height and width of the test section were not included in the 

foregoing because a detailed analysis, based on ideal vortices, has 

shown that the velocities induced by the bottom or sides were 
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negligible. Effects of surface tension on the instantaneous position of 
the vortex pair are deemed negligible and thus are not included as a 
parameter in equations (1) and (2). 

A dimensional analysis of equations (1) and (2) yields: 



is known as the Brunt-Vaisala Frequency. 

All parameters in equations (3) and (4) may be changed indepen- 
dently except r e /b 0 , which is taken as nature provides it. The primary 
reason for this is that a century of theoretical and experimental aero- 
dynamics research has been incapable of describing the details of the 
structure of the tip vortex to be used as the initial conditions in the 
viscous solution. It is surprising, but true, that until recently the 
importance of the generating surface shape and its influence upon 
both the initial tangential velocity profile and the initial turbulence in 
the vortex core had not been fully appreciated. Here the said influ- 
ence has been characterized in terms of an effective core radius with 
full awareness of its shortcomings. 




and 




in which 




(5) 
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B. GOVERNING DIFFERENTIAL EQUATIONS AND BOUNDARY 

CONDITIONS 

The generation of internal waves and the rise and demise of a 
vortex pair in a stratified medium may be analyzed through the use of 
the equations of motion for an incompressible fluid. These equations 
may be applied for both laminar and turbulent motions provided that a 
suitable turbulence closure model is used and the usual Boussinesq 
approximation (gravitational acceleration is much larger than fluid 
acceleration) is adopted. For the type of motion considered herein, 
the Boussinesq approximation is quite valid and has been used in the 
investigation of all types of internal waves in stratified fluids. 

For a two-dimensional flow, with y vertical and x horizontal, the 
Navier-Stokes equations of motion are: 



8u du 3u 1 3P - 0 

3t + u S +v Sy-=-p5T +vV2 u 



( 6 ) 



3v 3v 3v 1 3P _ 9 

3t + u K +v 3y = g'p37 +vV ^ 



(7) 



Differentiating equations (6) and (7) with respect to y and x respec- 
tively yields 

a 2 u a2 u a2 u a v au a 2 u i ap ap 1 a 2 P a ,_ 2 . 

3y3t + 3x3y + u 3x3y + 3y 3y + v 3 2 y ~ p^ 3y 3x p 3x3y + v 3y * * ' 



a 2 v a u av a 2 v a 2 v a 2 v 1 apap 1 a 2 P a ,_ 2 . 
a x at + ax ax + u ax 2 + axay + v axay ~ p^ ax ay p axay + v ax ( ^ 

Subtracting equation (9) from equation (8) yields 
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Subtracting equation (9) from equation (8) yields 



d fdu dv\ d [~ fdu dv\ 
dt [dy dx) + 3x dx J_ 



a 

+ 9y 



^av> 

ay ax. 




1 . 
~ ay ax 



i apap _ 2 fau av> 

F¥¥ + vV ¥"¥ 



( 10 ) 



i ap 



Solving equation (6) for — ^ yields 



I 

p dx 



= W’u-5| 



( 11 ) 



where 

d(].L) a a( j d( ) 

Dt ~ at u dx v ay 

1 dP 

Solving equation (7) for — yields 



1 aP Dv 

p jy = g + vV ^-Dt 



( 12 ) 



Substituting the results of equations (11) and (12) into equation (10) 
and defining vorticity as 



results in the following expression: 




( 14 ) 



When the gravitational acceleration is several orders of magnitude 
larger than the fluid accelerations, the Boussinesq approximation is 
appropriately invoked and the terms in brackets in equation (14) may 
be neglected. 

In order to deal with the case of a density stratified medium, the 
density is defined as 



where p 0 is the reference density, p (y) is the initial variation of density 
with y, and p’(x,y,t) is the fluctuating part of the density with time. 
Then 

§£. = ^£l 

3x dx 

and equation (14), neglecting the bracketed terms, becomes 



P = Po + p(y) + p'(x,y,t) 



- (15) 




(16) 
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The last term above represents the effect of the density gradient and 
gives rise to oppositely signed vorticity in a nonhomogeneous fluid. 
The diffusion of density is given by 



It + u fx + v 3y = vV2p 



(17) 



Substituting equation (15) for density in the above equation yields 



9p' 3p’ (dp 9p'^ (d 2 p' ^p 9 2 p’ > ) 

^ +u ^ + v (^ + ^J = v t&F + # + #J 



(18) 



The equation of continuity is 



9a 9v . 
dx + dy ~ 



(19) 



Adding equation (19) to the left side of equation (18) and simplifying 
results in 






3y 



ay 



9y 5 



( 20 ) 



Equations (16) and (20) are thus the governing differential equa- 
tions for the motion of a vortex pair in a density stratified fluid. 

For the development of boundary conditions, it is assumed at the 
outset that the influence of the sidewalls and bottom of the test 
section on the rise and demise of the vortex pair is negligible. The 
validity of this assumption, when used in conjunction with numerical 
computations, has been demonstrated by previous investigators [Refs. 
18-20]. The fluid domain can thus be considered to be bounded only 
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by a free surface which can be described as a function of x and time as 
follows: 



y = T|(x,t) 
where 



r|(x,o) = 0 



describes the initial undisturbed location of the free surface. 

The free surface requires both a kinematic and a dynamic bound- 
ary condition as described by Sarpkaya and Issacson [Ref. 22]. The 
kinematic condition states that any particle which lies on the free 
surface at any instant will never leave it. This leads to 



Dr| dT] 2r| 

Dt “ at + u ax “ v 



at y = r[ 



( 21 ) 



The- dynamic free surface condition requires that the pressure 
difference across the interface results in a force normal to the bound- 
ary which is due wholly to surface tension. This condition takes the 
form 



P = 



P a + a 




jn 

r 2 J 



( 22 ) 



where a is the surface tension, 4” and are the radii of curvature of 

Kl K2 



the free surface in any two orthogonal directions, and P - P a is the 
pressure difference across the interface. 
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As noted previously, the effects of surface tension can be assumed 
to be negligible. This observation has been verified experimentally by 
Gray [Ref. 17]. Thus, if the flow is considered inviscid, the pressure P 
within the fluid can be described by the unsteady Bernoulli equation as 
follows: 

90 a 2 P 

+ 2 + ^ (23) 



where 0 is the velocity potential such that 




v 



90 
" 9y 



(24) 



and q 2 = u 2 + v 2 . F(t) is an arbitrary function of time only. 

When the pressure just outside the liquid is constant (i.e., atmo- 
spheric), the free surface condition reduces to 



90 q 2 

af+2 + ^ = 0 at y = T1 



where F(t) has been included as part of 
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C. DIMENSIONLESS PARAMETERS 

It is convenient at this point to cast the governing equations, (15) 
and (20), in nondimensional forms, scaling each variable by a quantity 
characteristic of its expected magnitude. Two possible time scales 
exist. The dynamic time scale utilizes the time a characteristic length 
would be traversed by a fluid particle traveling at the characteristic 
velocity. The buoyant time scale is based on the natural buoyancy 
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frequency of the stratified flow, i.e., the Brunt-Vaisala frequency N 0 
defined by equation (5). Each of these scales gives a slightly different 
form of the normalized governing equations. 

1. Dynamic Scale 

Introducing U c and L c as the characteristic velocity and 
length, one has the following nondimensionalized quantities: 



Cm = CLc/Uc t m = Uct/Lc 



u m = u/U c 



v m = v/U c 



x m = x/Lc 



ym = y/Lc 



Pm - p/ Po 



(26) 



Substituting the above into equation (16) yields 



U c a 



Lc 5tnT ( ? 1 lc °) + L 3x m t UmUc lc 



a ( tt CmU c 



■) + iblb(; 



VmUc 



UU< 



4 vV “(Tr) + L^5^ lpop,nJ 



(27) 



Simplifying 

Uc 2 a^m Up 2 a(u m Cm) Uc 2 a(vmCm) Uc \7 2 r & dp'm 
Lc 2 atm Lc 2 dx m Lq’ dy m Lc a ^ 111 m Lc axm 



which becomes 



a Cm d(limCm) a(vmCm) V „ q r gLc ap' m 

atm + ax m + ay m U c Lc m + U? dx m 



or 



aCm a(u m C m) d(vmC m) _ J _„2 r 1 dp' m 
atm + 3x m + aym Re m m + Fv 2 5x m 



(28) 



( 29 ) 
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where 



Re = U c Lc/v (characteristic Reynolds number) 



F v = Us/VgLc (characteristic Froude number) 



Similarly, the density diffusion equation may be expressed in 
nondimensional terms by making equivalent substitutions into equa- 
tion (20) as follows: 




3(up') 
8x + 



3(vp‘) 

5y 




This becomes 



at m (p0pm) + Lc dx m (U c u mPoP'm) + (vmU c poPm) 

Uc v m 5, ,v \ — n . , ^ V (PoP'm) 

Lc ay m (popm) + L* V m (PoPnJ + L?— ^ 



Which upon simplification reduces to 

UcPo ap' m U c Po a(Ump'm) UcPo a(VmPm) 

Lc atm Lc axm Lc aym 



UcPo 

Lc 



v m 



apm vpo 
tym. Lj£ 




Pm + 



Vpo a 2 p m 

~W 






Multiplying both sides by Lc/poU c yields 



ap'm a(u m p'm) a(VmP'm) a pm V 

atm + ax m + ay m Vm ay m + UcLc 



Pm 



a2p 

ay. 



m 



( 30 ) 
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Incorporating the characteristic Reynolds number, this equation 
becomes 



dp'm d( u mP'm) d( v mP'm) 



at 



m 



dx. 



m 



dyi 



m 



_ d Pm 1 
m dym Re 



v 2 

V m p m + -^T- 



(31) 



The Brunt-Vaisala frequency given by equation (5) as 



- -'->1/2 



N 0 = 



z£ 32. 

Po dy 



may be written as 



N 2 _ _ fL 3( P°P m) . 

0 PoLc 3ym 

which becomes 

Kh 9 Pm 
g 3ym 

To further simplify equation (31), the following definitions are made: 



N* = g/Lc 



and 



(characteristic Brunt-Vaisala frequency 
squared) 



N| N|Lc 3p m 

n ” g - 3ym 



(33) 



(34) 



Substituting equation (34) into equation (31), the nondimensional 
form of the density diffusion equation in final form becomes 
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(35) 



dpm d(UmPin) ^( v mPm) 



at 



m 



ax 



m 



ayi 



m 



= n2 vm + r; 



Pm + 






— \ 
m 
2” 
m 



2. Buoyant Scale 

Again introducing U c and Lc as the characteristic velocity and 
length, the following nondimensional quantities are defined: 



u m 






Pm 



a JL 

Po 



In this case, however, time is nondimensionalized using the square of 
the characteristic Brunt-Vaisala frequency. Namely 



N? = g/Lc 



and then 

tm ~ tN c 
also 

Cm = C/N c 

Substituting the above into the governing equation given by equation 
(16) yields 

Nc ^7 + L^al^ ( u mU c CmN c ) + ( v mU c CmN c ) = 

(UNc ) +E ^--|-( PopW 
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simplifying 



rj o au , UcNc a 
c 9t m Lc 3x m 



(UmCm) + 



UcNc_3_ f , ; _ N c _ 2 . 

Lc ^ym (Vm ^ m) " IJ vV mCm + 



g dp’m 
Lc ax m 



which becomes 



dCm Uc ^(UmCm) Uc a(VmCm) V \j2r & ^P'm 
3t m + LcN c ax m + L C N C dy m ~ N^L? m ^ m + ax m 



136) 



but 



n|^ - n|" 1 

Thus equation (36) becomes 



a Cm Uc fd(UmCm) aCVmCm)' ) V Uc jn2r ^Pm 

5tm + -x^gLc v ^ Xm + dym J ^[gL^L c Uc m m + ^ x m 



or 



aCm ^ ra(umCm) a(v m ^ni)l F v „ 0y ap'm 
3t m V L 3x m + dy m J-Re V m^m+ 3Xm 



(37) 



where F v and Re are as defined above. 

The density diffusion equation is similarly nondimensional- 



ized as 



ap' a(up') a(vp') 
at + ax + ay " ~ 



v 
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which becomes 



9 1 9 19 

Nc ( P°P ,m) + (umUcpop'm) + (vmUcPoPm) 

U c Vm 9,_ x \ „ 0 , ,, v 9 2 , _ * 

" Lc 3y m (poPm) + I? V m (PoPnJ +IY^(PoPm) 



The above equation can be simplified to yield 

» T dp' m . U c po 9 , % . UcPo 3 , , ^ 

PoNc atm Lc 9x m (u mPm) + ^ ^Mm) 

~UcPo _ ap m p 0 V „ 9 Vp 0 a 2 pm 

T Vm ay^ + v ” pm + Ti"5S" 



Introducing F v and Re as before and recalling that N 2 = one has 



dp'™ , F 
at, +Fv 



Lm 



a(UmP’m) 


a(v m p'm)l 


. dx m 


ay m J 



•c' ^Pm 

FvVm ay m + 



FV 

Re 



V m p m + 32 






'm 



(38) 



Now as developed above in equation (34) 



n 2 = 



ap 



m 



where n 2 = 



NJ N|Lc 

w~ 



3ym “NJ- g 

a p m 

Substituting for in equation (38) yields the nondimensionalized 
density diffusion equation 

ap'm 



at 



+ Fv 



m 



ra(UmPm) 


a(v m Prn)l 


L ax m 


J 



,, o Fv 
F v n^v m + ^ 



„ o _ 1 a 2 p m N ' 

V mP m + - P" 

y m J 



(39) 
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Equations (37) and (39) are valid when F v «l and buoyancy 
dominates the flow. When F v approaches zero, the equations that 
result from equations (37) and (39) describe the propagation of linear 
internal waves. The buoyant scaled forms of the governing equations 
are of interest to the present investigation because, for submerged 
bodies of naval interest, F v is typically about 0.001 and thus signifi- 
cantly less than one. 

Having established the buoyant scale as the form of interest, 
the boundary conditions can also be nondimensionalized using the 
same technique. From equation (21), 



ar +u s? =v at y =, i 



+ u- 



which becomes 




w ' mi ~ v, — Hi ~ ■ 1 111 — - 

9tm "** Lc 9xni c m 



Uc u m 



where 



Tim = Tt/L c 



Simplifying the above yields 



dTlm Uc 9rim Uc 



9t m N c Lc Um 9xm N^Lc m 



Introducing F v as before 



9t m + Fv Um 9x m " FvVm at Ym = 1,111 



(40) 
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and 



0 

0 m = (nondimensional velocity potential) 



Equation (25) can be nondimensionalized as follows: 

,2 



30 Q~_ „ 

•sr + ^+gn = 0 



becomes 



N< 



3(0 m N c L?) U?q2 



at 



ra 



+ 2 + “ 0 



where 



<lm = + v m 



Further simplifying equation (41) yields 
30m ~ U? q^j gT|m « 

■5t^- + N^— + N^=° 



but 



N^Lc 



= 1 



and 



U? „ 2 
N^i= F ? 



Thus equation (42) becomes 
30m 






3t 



m 



+ Fy 2 + Tim - 0 at ym - Tim 



(41) 



(42) 
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For the purposes of the present investigation, the flow will be 
assumed to be inviscid, i.e., the viscous diffusion will be ignored to a 
first order approximation. The nondimensionalized forms of the 
governing equations and boundary conditions are summarized in 
Figure 1. 
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HI. NUMERICAL SOLUTION TECHNIQUES 



A. INTRODUCTION 

The significant difference between this study and the work of 
previous investigators is the introduction of a deformable free surface 
and the resulting complex interplay between kinetic and potential 
energies. The computational domain thus has an unknown boundary 
on which a double condition must be imposed as previously outlined. 
This complexity, as well as several other specific features of this 
problem, directly influences the ability of a numerical method to con- 
verge to an adequate solution. 

The governing differential equations and boundary conditions are 
both nonlinear. Sarpkaya and Issacson [Ref. 23] note that for small 
amplitude waves (amplitude « wavelength), the boundary conditions 
at the free surface may be linearized. This approach, although possibly 
valid in the vicinity of surface striations, would not be valid for the scar 
front where observed light diffraction patterns [Ref. 17] indicate very 
small radii of curvature. Additionally, Haussling and Coleman [Ref. 24] 
have demonstrated numerically the importance of nonlinear terms in 
solving the free surface problem in the vicinity of the generation 
source. This is the case encountered in the region of the scar front. 

The rate of change of the free surface deformations also tends to 
be slow. The scar pattern develops as the vortex pair rises to its 
maximum height and then the scar is trapped by and slaved to the 
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vortex beneath the surface. The scars and striations are thus not small 
amplitude surface waves propagating independently from their 
generator, but are in fact local disturbances whose generation, growth, 
and propagation are controlled by the vortices. This aspect, observed 
experimentally by Gray [Ref. 17], separates this problem from other 
research in the field of deformable free surfaces. Previous investiga- 
tors, applying numerical methods to such surfaces, have dealt almost 
exclusively with surface waves generated by a moving body (surfaced or 
submerged) and propagating away independently. 

The numerical methods reviewed for solution of this problem can 
be broadly categorized into Finite Differences, Finite Elements, 
Boundary Integral Equations, and Hybrid Methods which employ com- 
binations of the others. For the purposes of this study, such catego- 
rization is based on the manner in which the governing equations are 
tackled in the computational domain. It is noted that all of the 
methods' reviewed utilize some form of finite difference scheme with 
respect to time. The applicability of each method and the associated 
advantages and disadvantages are discussed and summarized at the 
end of this chapter. 

B. FINITE DIFFERENCE TECHNIQUES 

The use of Finite Difference techniques in solving partial differen- 
tial equations has been well established and successfully implemented 
by several investigators. For the case of a nondeforming free surface, 
the governing equations are solved in the computational domain using 
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dju g + Mi. vV 2, 

3t + 8x dy ^ Po dx 



( 30 ) 



and 



3p* 8(up) 3(vp) 

3t + 3x + dy 



- V f- + wv + v|| 






and the Biot-Savart integrals 



u(x,y) = J 



(31) 




(x -X’) c (x’, y’) dx'dy* 
2xr 2 



(32) 



where 

r 2 = (x - x’) + y - y’) 

x\ y': position of the vorticity 

x, y: point where u and v are to be calculated 

In nondimensional form, the buoyantly scaled forms of equations (31) 
and (32) become 
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u m ( x m>ym) — 



* (y’m "Yin) Cm (x'm. y'm) dx’ m dy’m LcN c 

2^ U c 

J 

A 



v m (Xm.ym) = 



' (x m -x' m ) Cm (x'm. y'm) dx’ m dy’ m LcN c 

2^ U c 

J 

A 



Upon simplification, the equations become 



Fy Um (Xm.ym) - 



(y'm “ ym) Cm (x'm. y'm) dx' m dy m 



27tr m 



(33) 



Fy v m (Xm.ym) — 



(x m -x' m ) Cm (x’m. y'm) dx' m dy’ : 



m 



2xr 



m 



where parameters are nondimensionalized as before. 

The' governing equations can then also be rewritten as 

3 Cm 3(F V Um Cm) 3(F V v m Cm) F v - 3p'm 

3t m 3x m 3y m Re m m 3x m 



(34) 



and 



3p'm 3(F V U m P m) 

3t m + 3x m 



3(F v v m p'm) 



3y : 



m 



o Fy 
Fyv m n 2 + ^ 



- 



X 72 > , ^ 2 Pni 

pm 



where F v u m and F v v m are calculated directly using the nondimension- 
alized forms of the Biot-Savart Integrals equations (33) and (34). 
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The treatment of the governing differential equations in a stream 
function formulation is thus relatively straightforward. The complica- 
tion is associated with the determination of the free surface location in 
conjunction with the solution of the governing equations, all preferably 
simultaneously, as noted by Yeung [Ref. 25]. 

Finite Difference methods are most suitable, or at least simplest 
to implement, for a boundary geometry that is rectilinear. The 
deformable free surface requires the use of “irregular stars" for differ- 
ence schemes near the fluid boundary. These “irregular stars” are of 
an inferior accuracy when compared to the remainder of the grid. 
Thus, either the mesh must be refined or the difference formula 
changed to a higher order if accuracy is to be maintained near the free 
surface and especially in the vicinity of the scars. Yeung [Ref. 25] 
notes that, with the location of the free surface unknown a priori, we 
have the unfortunate situation that the regions that demand the great- 
est accuracy are precisely those where it is hardest to achieve. Also, 
since the shape of the free surface affects the migration of the vortex 
pair, a loss of accuracy in determining the free surface will be 
reflected in the determination of the location of the vortex pair. This 
complication is one not included in the investigations reviewed below 
where the generating body moves independently of the free surface 
location. 

Conformal transformations can simplify one aspect of the problem 
by transforming the real geometry with a deformed free surface to a 
rectilinear mesh. The boundary conditions are thus simplified but the 
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resulting field equations are increased in complexity. Haussling and 
Coleman have successfully utilized transformation functions of this 
form to generate a time independent computational region in which 
nonlinear free surface boundary conditions are applied. 

The Finite Difference method primarily benefits from the direct 
solution, of the governing differential equations as opposed to- the 
development of intermediate integral relations required in the inte- 
gral boundary techniques. The cost of such an approach is in the 
computational intensity required to iteratively solve the finite differ- 
ence forms over the computational domain. Haussling and Coleman 
have demonstrated the use of a successive over relaxation technique to 
solve steady state nonlinear free surface boundary condition problems 
as discussed above. 

A modification of successive over relaxation, and one offering 
increased rates of convergence, has been demonstrated by Brandt, 
Dendy, mid Ruppel [Ref. 26]. Their technique utilizes a multigrid 
solution which solves for low-frequency components of error on a 
course grid where the calculation is relatively inexpensive, and high- 
frequency components of error on a fine grid where successive over- 
relaxation is efficient. 

Theodussiou and Sousa [Ref. 27] also utilized a modified grid sys- 
tem to speed convergence by staggering their grid structure such that 
pressure was defined at the center of the discretized control domain 
and velocity components were defined at the center of the control 
domain faces. This arrangement has the convenient feature that the 
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velocity components are stored at just the points at which they are 
required for the calculation of their advective contributions. The 
pressure gradients can be represented by central differences without 
inducing non-physical oscillations in the pressure distribution. Note, 
however, that both of these multigrid techniques add to the complex- 
ity of establishing a mesh which moves with the deforming free -sur- 
face. Both sets of authors suggest, however, that this added complex- 
ity is outweighed by the savings in ease of convergence. 

Ohring [Ref. 28] and Ohring and Telste [Ref. 29] have also partially 
circumvented the computational intensity of an iterative technique by 
directly solving the finite difference equations resulting from the 
Laplace Equation. The technique employed utilizes a fourth-order 
solver to diagonally decompose the resulting coefficient matrix. Taylor 
series approximations are also suggested for application to the nonlin- 
ear free surface boundary conditions. It is suspected, however, that 
the computational benefits derived from such a technique will be lost 
when the coefficient matrix becomes nonlinear as would result from 
the governing equations in this problem. 

For completeness, as noted previously, finite differencing in time 
is required for all the numerical methods reviewed. Sarpkaya and his 
students [Refs. 17-20] have successfully employed an upwind -differ- 
encing scheme in time and verified it with experimental results. 
Yeung notes that the free surface conditions, being first order in time, 
can be used to advance the solution of the elevation and velocity 
potential on the free boundary. However, the difference form utilized 



40 



can dramatically affect the stability and thus a modified Euler method 
is suggested since it is unconditionally stable. Haussling and Coleman 
successfully utilized this technique for advancing their solution in 
time. 

C. FINITE ELEMENT METHODS 

The Finite Element and Finite Difference methods share the 
common feature that both attempt to solve the governing differential 
equations directly, differing only in methodology. The Finite Element 
Method is one based on the method of weighted residuals. The usual 
procedure consists of first subdividing the domain of interest into a 
mesh of finite-sized subregions, within each of which the solution is 
represented by some convenient choice of trial functions, usually 
polynomials. The trial functions are determined by substituting them 
into the governing equations and requiring the integrated error or 
residual based on certain weighing functions to vanish. An integration 
by parts is normally preformed to reduce the interelement continuity 
requirements of the trial functions and to incorporate the 
nonhomogeneous boundary conditions. The weighing functions, also 
known as test functions, can be chosen in a variety of ways. A “weak 
formulation,” such as the Galerkin Method, makes the space of the 
test functions identical to that of the trial functions. In contrast, a 
“strong formulation” is one based on the existence of a variational 
principle where a functional is made stationary. Specifics of finite 
element techniques may be found in Zienkiewicz [Ref. 30] and Dhalt 
and Touzot [Ref. 31], 
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The nonlinear aspects of this problem affect the Finite Element 
Method in much the same way as the Finite Difference Method. The 
coefficient matrix of a representation such as 
[K] [U] = [B] 
where 

[K] = coefficient matrix 
[U] = matrix of unknowns 
[B] = forcing function matrix 

is neither symmetric nor independent of [U], as in the case of a linear 
problem. An iterative technique is thus required to solve this problem 
at a given instance in time. As a result, much of the advantage of 
reduced storage and computation time inherent to finite element for- 
mulations is lost. The choice of iterative techniques does not differ 
from that available to Finite Difference Methods and, thus, once the 
computational domain is discretized, no significant difference in 
problem solution is involved. 

The discretization of the computational domain is an area where 
the Finite Element Method does have distinct advantages. The intro- 
duction of curvilinear or isoparametric elements of higher order 
allows one to cope with any arbitrary boundary geometry with little 
loss of accuracy. Yeung notes that, although this particular advantage 
is reduced when Finite Difference Methods are used with conformal 
transformations. Finite Element Methods still retain superiority in 
flexibility, particularly in the case where varying size and shape 
elements are introduced to overcome local irregularities. Such 
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irregularities are exactly the problem involved in the vicinity of the 
deformable free surface. Curvilinear elements could be used to fit the 
free surface with finer elements incorporated in the vicinity of the 
scar front. This could be accomplished with little loss of accuracy and 
less complexity than that which results from the use of “irregular 
stars" in the Finite Difference Method. Note, however, that, since the 
scar front moves with time, an adaptive mesh refinement technique 
would be required to appropriately place the finer elements in the 
vicinity of the scar. Sarpkaya and Hiriart [Ref. 32] successfully utilized 
varying size elements in the vicinity of the free surface in conjunction 
with their moving net computation. Although the free surface in their 
case differed from that involved in this problem, the basic concept of a 
“flexible mesh" remains unchanged. 

Admittedly, the basic problem of determining the location of the 
free surface in conjunction with solving the governing equations 
remains' one of trial and error, regardless of the type of discretization 
employed. Larock.and Taylor [Ref. 33] adjusted the free surface loca- 
tion to achieve tangency to the surface velocities calculated based on 
an assumed free surface position. The pressure boundary condition 
was then used as a check of this resulting location. Larock and Taylor 
note, however, that such a technique will not work well where high 
curvature or Froude numbers (F v ) less than one are involved, as is the 
case in the present investigation. As an alternative, Sarpkaya and 
Hiriart adjusted their free surface location to satisfy both the pressure 
and velocity boundary conditions simultaneously. 
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The Finite Element Method, when formulated using the varia- 
tional techniques, also benefits from the fact that the boundary condi- 
tions may be included as an integral part of the functional rather than 
dealt with separately. Bai and Yeung [Ref. 34] and others have suc- 
cessfully employed variational techniques to directly incorporate 
boundary conditions in the solution of linear water wave problems^ 

D. BOUNDARY INTEGRAL EQUATION METHODS 

The term “Boundary Integral Equation Method” can be applied to 
a large group of numerical techniques that include Green’s Function 
formulations. Spectral Methods, and a boundary element application of 
the Finite Element Method. In all cases, the solution approach differs 
dramatically from that of Finite Differences and Finite Elements in 
that the governing equations are not attacked directly. Instead, the 
problem is solved by satisfying boundary conditions on a discretized 
boundary and thus reducing the spatial dimensions of the problem by 
one. Physical quantities such as wave height and fluid pressures are 
required and solved for only on the boundaries. Interior data, although 
available based on the boundary solution, is not specifically required. 
Boundary Integral Equation Methods thus have the distinct advantage 
that only the physical quantities specifically desired on the free sur- 
face are required for the solution. 

Basically, there are two approaches to boundary integral formula- 
tions. Either an approximating function is chosen on the boundary 
that satisfies the governing equations in the domain and approximates 
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the boundary conditions, or vice versa. In both cases, the solution 
technique can be further divided into indirect and direct methods. 

In indirect approaches, the fundamental solution is approximated 
on the boundary by a function with unknown coefficients. These coef- 
ficients are found by satisfying the boundary conditions. Distributed 
singularity methods, such as simple-source distributions, are an exam- 
ple of this approach. In this case, the fundamental solution is repre- 
sented on a discretized boundary by distributed simple sources of 
unknown strength at known locations. The strength of each source is 
then determined based on the solution of the boundary conditions. 

Direct methods find the fundamental solution through the use of 
Green’s Function formulations, which directly incorporate the 
governing equations. The resulting Green’s Function integrals are 
typically solved using quadrature or point kernel techniques. The 
main disadvantage in this approach is that, for the governing equations 
involved“in this problem, there is no straightforward Green’s Function 
formulation that leads to fundamental solution. 

The main advantage of the Boundary Integral Equation Method is, 
then, the ability to directly discretize the free surface without a loss of 
accuracy. Information is thus obtained exactly where required. How- 
ever, Brebbia [Ref. 35] notes that this is not without the sacrifice of a 
symmetric coefficient matrix such as that which is common to the 
Finite Element Method. Yeung notes that, in general, the great 
reduction in the size of the matrix outweighs the added complexity of 
solving a nonsymmetric system of equations. 
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The disadvantages of Boundary Integral Equation Methods are 
directly associated with errors that are the result of discretization. 
This discretization plays such a large role in the formulation of an 
efficient boundary integral equation that Brebbia refers to these tech- 
niques as “Boundary Element Methods.” In all formulation tech- 
niques, -elements of one type or another are formed over which a 
fundamental solution must be approximated either by distributed sin- 
gularities, Green’s Function integrals, or trial functions for finite 
elements. 

Discretization errors result both from collocation errors and geo- 
metric surface errors. Collocation errors are primarily the result of 
leakage, as noted by Hunt [Ref. 36]. Since the boundary conditions are 
satisfied exactly only at discrete points, the remainder of the free sur- 
face is “porous” by comparison. Leakage errors are reduced by 
increasing the number of free surface elements and thus the number 
of discretization points. This, of course, is done at the expense of 
computational intensity. Geometric surface errors result from the 
approximation of a curved free surface by linear elements. This prob- 
lem becomes particularly noticeable in areas of high curvature, such as 
scar fronts. Higgins and Cokelet [Ref. 37] noted that this problem can 
be partially circumvented by using a Lagrangian description of 
“marked” particles on the free surface. These particles tended to 
concentrate in the regions of highest curvature, thus giving improved 
accuracy exactly where needed. 
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Finite Element Method formulations using boundary elements 
only can also circumvent some of the errors associated with the collo- 
cation and geometric surface by using curvilinear or higher order 
parametric one-dimensional elements. The problem then is one of 
finding suitable fundamental solutions that satisfy both the governing 
equations and continuity requirements between elements. 

Finally, Boundary Integral Equation Methods are complicated by 
the necessity to include nonlinear terms in the boundary condition for 
the free surface. As previously noted, this is necessary in order to 
maintain accuracy in the vicinity of the scar front where very small 
radii of curvature occur. Several researchers have successfully applied 
various forms of the Boundary Integral Equation Method to Laplace 
equations with linearized free surface conditions, but only a few have 
incorporated a nonlinear condition. Faltinsen [Ref. 38] incorporated 
nonlinear conditions by using the properties of the fluid particles on 
the free surface at one instance in time to establish a new free surface 
location and step the solution forward in time. 

E. HYBRID METHODS 

The use of a different technique in different portions of the com- 
putational domain results in a hybrid numerical formulation. This type 
of formulation attempts to maximize the benefits of any one particular 
method by employing it only in regions where its accuracy remains 
high. A secondary objective is to reduce the computational intensity of 
the overall routine by using a coarser, less time-intensive technique in 
areas where accuracy is not of particular concern. This is exactly the 
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case in the regions far removed from the free surface. The usual 
approach is to take advantage of the availability of analytical solutions 
in regions where the flow geometry is relatively simple. This 
approach allows the number of mesh points to be reduced with a 
resulting decrease in required storage and computational intensity. A 
further advantage is that the analytical solutions can be chosen to 
permit a simple solution of radiation type boundary conditions, as 
noted by Yeung [Ref. 25]. 

The disadvantages of Hybrid Methods are associated with the 
necessity to properly match the solutions of different subdomains in 
the overall computational domain. This matching may result in 
numerical perturbations to the solution technique if not properly 
employed. Additionally, the advantage of reduced total grid points is 
somewhat offset by the added computations required to match 
solutions at the common boundaries. 

Yeung notes that “the more successful Hybrid methods have so far 
been restricted to linearized problems where analytical solutions in 
the exterior regions could be obtained without too much difficulty. In 
particular, treatment of steady flows in a uniform stream or time-har- 
monic flows with linearized boundary conditions have been quite well 
established." [Ref. 25] 

Bai and Yeung [Ref. 34] utilized a finite element grid in the vicinity 
of the free surface and an eigenfunction solution in the exterior region 
to reduce the computational intensity of their routine. Chang and Pien 
[Ref. 39] used a source distribution to formulate a Boundary Integral 
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Equation Method near the free surface and a finite difference routine 
in the exterior regions to calculate the hydrodynamic forces on a body 
moving beneath a free surface. It should be noted, however, that both 
sets of investigators used a Laplace formulation with linearized bound- 
ary conditions to obtain their results. 

For ..completeness, it is noted that Hybrid methods could be-£on- 
sidered to include various numerical techniques that employ the same 
method throughout the computational domain but in different man- 
ners in each subdomain. Such is the case, for example, when a finite 
element method is employed with a varying grid and/or element type 
in various regions of the computational domain. 

F. NUMERICAL METHODS, CONCLUSIONS, AND SELECTION 

The synopsis of advantages and disadvantages for each numerical 
method listed in Table 1 provides a good source of information for 
making a wise choice of formulation to be used in this problem. The 
“best” method is one which will meet the goals of the investigation 
while maximizing the advantages in areas of particular concern. 

In this case, accuracy in the vicinity of the free surface is of par- 
ticular importance since ultimately it is the free surface shape which 
is desired. This accuracy must be obtained with full consideration of 
the necessity to include nonlinear boundary conditions while mini- 
mizing the computational intensity of an iterative process. Tuck, in a 
paper by Bai and Yeung, observed that “to a certain extent the ‘best* 
method will always be that which appeals most to the person pro- 
gramming it, and hence that which gives him the greatest chance of 
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writing a successful program irrespective of efficiency. There is no 
less efficient program than one which does not work at all. n [Ref. 34] 
Given that a successful program can be developed (a pretty big given), 
it is then necessary to maximize both efficiency and accuracy. 

With full consideration of the concerns given above, the Boundary 
Integral- Equation Method surfaces as the method of choice for the 
following reasons: 

1. Computational intensity is minimized by reducing the spatial 
dimensions by one. This aspect is particularly noticeable when 
consideration is given to the iterative requirements of this 
problem. This substantial improvement directly incorporates 
any similar advantage that can be gained by using a varying Finite 
Element Method or Hybrid Method. 

2. Accuracy is maintained at the free surface by incorporating the 
best aspects of Finite Element Methods and a Lagrangian 
description of marked particles. This direct discretization of 
the free surface provides both accuracy where needed most and 
an exact solution of boundary conditions at the marked points. 

3. Storage requirements are drastically reduced since only physical 
quantities at the free surface are required or needed. 

The implementation of the Boundary Integral Equation Method 
through the use of distributed vortices will be outlined in detail in the 
following section. 
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COMPARISON OF NUMERICAL METHODS 
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erning dmerential equations. tionally intensive when 

. The boundary conditions can be employed with an iterative pro- 

included as an integral part of cedure (nonlinear coefficient 

the formulation when employed matrix. , 

with variational techniques. 
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cable, can simplify radiation- rnay be reduced when employed 

type boundary conditions. in conjunction with other for- 



IV. DISTRIBUTED VORTEX MODEL 



A. INTRODUCTION 

Discrete vortices, with or without a core, or vortex sheets have 
been used as boundary elements to simulate separated and unsepa- 
rated flows. The method consists of the determination of the appro- 
priate strength of the vortices at each time interval through the use of 
the governing equations. The vortices are then convected to their 
new positions and the process is repeated. In spite of its simplicity, 
the distributed vortex model presents several difficulties, all of which 
are related to discretization and the use of vortices. The evaluation of 
the governing integral equation cannot accurately be accomplished 
merely by applying a standard integration formula. A vortex sheet or a 
string of point vortices is unstable to small sinusoidal disturbances of 
any wavelength. This phenomenon of Helmholtz instability persists in 
curved nonuniform vortex sheets, at least for short waves, unless the 
sheet is rapidly stretching. In other words, the round-off and trunca- 
tion errors are rapidly amplified to cause the chaotic motion which 
often ruins practical calculation. If it is granted that it is the growth of 
short waves which can ruin calculations with vortex sheets, it is sensi- 
ble to consider ways of removing the instability. This is because the 
instability is introduced by the step of replacing a shear layer of small, 
but finite, thickness by a vortex sheet. One could give up the vortex 
sheet approximation and return to the computation of the evolution of 
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a thin layer. This is without doubt the most satisfactory procedure, but 
it involves much more computation. 

An alternative approach is to modify the governing differential 
equation to allow for finite thickness but the resulting equation, while 
only a little more complicated than the original equation, has not 
proved amenable to computation. 

Another possibility is to apply a linear smoothing formula, such as 
that introduced by Longuet-Higgins and Cokelet [Ref. 37] in their 
work on nonlinear water waves. The subsequent applications, includ- 
ing the one discussed herein, have shown that chaotic motion sets in 
sooner or later regardless of the smoothing. The repositioning tech- 
nique introduced by Fink and Soh [Ref. 40] and used subsequently by 
Sarpkaya and Shoaff [Ref. 41] removes the most unstable mode and 
reduces the growth of the higher modes of Helmholtz instability. 
However, it does not prevent the growth of spurious waves along the 
vortex sheet. The disadvantage of the smoothing, either through the 
use of a numerical filter or through the repositioning of the vortices, is 
that it is not clear in general what the relationship is between the 
results achieved and the unknown exact solution. The problem may- 
not possess a solution for all time, and in this case the use of smooth- 
ing could yield an acceptable-looking solution where none in fact 
exists. Alternatively, the solution arrived at through smoothing may 
not be even close to the exact solution. 
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B. APPLICATION OF THE DISTRIBUTED VORTEX METHOD 



The liquid surface was represented by a number of point vortices. 
Their positions were symmetrical with respect to the x axis, situated 
on the free surface. However, the strength of a vortex at (x, y, nor- 
malized by b Q ) was opposite to that of a vortex (image vortex) at (-x, y) . 
From a mathematical point of view, one would like to have the vortices 
extend from -~> to +<» and the two trailing vortices originate at (1/2, 
-°o) and at -1/2, -®o). This is impossible from a numerical point of 
view. Observations have shown that the free surface rises in a rela- 
tively small region directly above the trailing vortices. The remainder 
of the free surface remains undisturbed. These observations and 
several sample calculations led to the conclusion that the free surface 
can be restricted to a region extending from x = -10 to x = +10. 
Furthermore, the trailing vortices are assumed to originate at a depth 
of y = -5. 

The complex function representing the two trailing vortices 
below the free surface and the vortices on the free surface is given by 



irm 

2tc 



In (z + zi) 



in which the strengths of the free-surface vortices are normalized by 
the strength of the trailing vortex. All coordinates are normalized by 
b 0 - Then the velocities, normalized by V 0 = (r 0 /27tb 0 ) anywhere in the 
fluid medium is given by 
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( 35 ) 



dw ipo 1 ipo 1 iFm 1 iTm 1 

U - IV = — = - TTT — Z - “ + 7TZ" Z - + 



dt 2k ( z _z 0 ) 2k ( z + z 0 ) 2k (z- z t ) 2 k ( z+ Zi) 



The normalized boundary condition at the free surface is given by 
2 



D0 m Qm T lm 

v 



Dt m 2 + W ~ 0 



(36) 



Thus, at least theoretically speaking, one can calculate the poten- 
tial function 0 from the complex velocity potential w (the real part of 
w), the velocity of the vortices from equation (35), the elevation of the 
free surface from equation (36) for a given Froude number Fv (Fv = 
V 0 /Vgb 0 ), and trace the evolution of the free surface as a function of 
time. This relatively simple-sounding procedure is anything but sim- 
ple, primarily because of the fact that the numerical instabilities even- 
tually lead to large-scale instabilities, as noted in the introduction to 
this section. The use of various smoothing techniques was a viable 
option and, in fact, was tried at various stages of the calculations. The 
results have shown that the growth of the instabilities increases with 
decreasing Froude number. Neither the use of smoothing techniques 
(e.g., the Longuet-Higgins technique) nor the use of vortex sheets, 
vice discrete vortices, can eliminate the eventual development of a 
chaotic behavior on the free surface. It is because of this reason that 
the use of smoothing was ruled out and the calculations were 
restricted to a relatively high Froude number (Fv = 1.125). 

The specific details of the numerical calculations are as follows: 
(1) assign the position of the vortices; (2) find the strength of the vor- 
tices and the velocity potential assuming the free surface to be rigid; 
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(3) advance the position of the trailing vortices for a time interval t 
(0.01 in the calculations); (4) recalculate the velocity of the surface 
vortices and advance them forward in time for a time interval t; (5) 
calculate the strength of the free surface vortices by iteration until the 
free surface condition is satisfied; (6) calculate the velocity of the sur- 
face and trailing vortices and advance their positions for a time inter- 
val t; and (7) repeat the calculations by returning to step (5). 

The procedure described above is relatively simple and does not 
require excessive computer time (about 20 minutes on IBM 3033). 
However, the free surface eventually does become chaotic. Some of 
the instabilities can be alleviated without smoothing through a judi- 
cious selection of the initial position of the surface vortices. Numerous 
calculations have shown that the vortices near the y axis begin to slide 
sideways as the free surface (or the vortex sheet) stretches. Conse- 
quently, the thinly populated regions of the sheet do not yield a 
smooth free surface and the flow begins to leak between the vortices. 
Also, the regions where the free surface is depressed (where the scars 
develop) become overpopulated with vortices, leading to the growth of 
short-wavelength Helmholtz instability. This problem can easily be 
alleviated by packing the vortices more densely near the y axis at the 
start of the calculations so that the entire surface becomes more or 
less uniformly represented at later times. This simple technique has 
been used in the results presented herein. An exponential function 
was employed to assign the initial vortex spacings so that near the 
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y axis the vortices were separated by a distance of about 0.03 and by a 
distance of about 0.25 towards the end of the vortex sheet. 

C. RESULTS OF THE NUMERICAL CALCULATIONS 

The evolution of the free surface with time is shown in Figures 2. 1 
through 2.27 for Fv = 1.125. The normalized time in these plots is 
given by T = (V Q t/b 0 - d 0 /b 0 ) where d 0 is the initial depth of the 
trailing vortices. The time T = 0 corresponds to that at which the 
trailing vortices would have reached the undisturbed free surface had 
they continued to move with their initial mutual induction velocity. 
The use of other normalized times is not suitable since they tend to 
depend on the initial position of the trailing vortices. 

Figures 2.1 through 2.27 show that the free surface directly above 
the vortices rises rapidly while the adjacent portions of the surface 
depress and form two strong scars. Unlike the rigid surface case, 
where the trailing vortices continue to move towards right and left at a 
depth of about y = -0.5, the trailing vortices below the deforming sur- 
face almost come to rest at a point slightly above the free surface. 
Furthermore, their initial spacing remains nearly constant. This sug- 
gests that the Kelvin oval formed by the trailing vortices pushes the 
free surface up as if it were rigid during its upward migration. Evi- 
dently, this finding is based on the assumption that the trailing vor- 
tices are point vortices and are not subjected to viscous and turbulent 
diffusion. In reality, the vortices quickly become turbulent, their core 
radius increases, and the vorticity is diffused over an ever-increasing 
area with the passage of time. The amount of diffusion depends on the 
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initial position of the trailing vortices. Consequently, the trailing vor- 
tices emanating from large depths diffuse over a larger area relative to 
those emanating from smaller depths. There is at present no suitable 
mathematical means to deal with the turbulent diffusion of such vor- 
tices. The best one can do is to generate the trailing vortices at 
depths sufficiently close to the free surface to prevent excessive dif- 
fusion and yet sufficiently far so that the free surface remains 
undeformed at the start of the motion. Extensive calculations and 
experimental observations have shown that the free surface does not 
begin to deform until the trailing vortices reach a depth of about y = 
-1. Thus, it is perfectly safe to place the trailing vortices at y = -3 at 
the start of the calculations. 
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V. EXPERIMENTS 



A. EXPERIMENTAL APPARATUS AND PROCEDURES 

Experiments were conducted in a water basin. It consisted of a 
12-foot-long, 3-foot-wide tank with aluminum bottom and plexiglass 
walls (see Figure 3). Additional equipment consisted of plumbing for 
the filling and emptying of the tank, a collimated light source, and the 
dye system for flow visualization. 

The two-dimensional trailing vortices were generated through 
the use of two counter-rotating plates (see Figure 3). The mechanism 
rotating the plates was located at the bottom of the tank and below a 
plexiglass plate spanning the entire tank. In other words, the motion 
of the driving system did not disturb the flow above the plexiglass. 
The mechanism was actuated by releasing a weight attached to the 
common driving shaft. Fluorescent dye was introduced slowly into the 
region between the plates through the use of two holes connected to 
two dye reservoirs. The reason for the use of two reservoirs was that 
different colors of dye can be used to visualize each trailing vortex. 

B. PROCEDURES 

Experiments were initiated by filling the tank to a suitable level, 
removing any air bubbles, bringing the plates to their full open posi- 
tions, waiting for a sufficient period of time for the fluid to come to 
rest, introducing the neutrally buoyant dye, actuating the light source 
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and the video system, and initiating the rotation of the plates by sud- 
denly releasing the load attached to the driving shaft. The plates 
rotated smoothly and then came to rest on the plexiglass. In other 
words, the plates generating the vortices literally “disappeared.” 
Thus, no other vortices were generated. It is a well-known fact that 
this is not the case with piston-generated vortices. When the piston 
stops, two additional vortices (from a two-dimensional piston) or 
another vortex ring (from an axisymmetric piston) are generated. The 
mechanism used in the present investigation effectively prevents the 
generation of secondary vortices by simply disappearing. 

The trailing vortices rise under their mutual induction velocity, 
quickly give rise to a Kelvin oval, and smoothly approach the free sur- 
face. When the vortices are at a distance of about 1.0 from the free 
surface, the free surface begins to rise and forms a smooth hump, bor- 
dered by two scars. For the Froude numbers achievable (maximum 
0.35), the trailing vortices turn quickly to right and left (as if they 
were approaching a rigid surface), the scar front moves in the respec- 
tive directions ahead of the vortex, and the hump between the vortices 
begins to recede. During the later stages of the motion, the scars are 
slaved to the vortices and continue to move ahead and in the direction 
of the vortices. 

Figures 4.1 through 4.12 show the evolution of the free surface at 
various times T for a Froude number of Fv = 0.35. The vortex centers 
are clearly visible in these pictures because of the additional assistance 
provided by the dissolved air bubbles to the flow visualization efforts. 
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VI. DISCUSSION OF RESULTS 



The computer code based on the boundary element methods 
through the use of the discrete vortices has been used to carry out 
calculations for various Froude numbers in order to predict the evolu- 
tion of the free surface deformation. The type of the vortex distribu- 
tion, time increment, and the iteration techniques have been varied 
within limits to explore the differences in the motion of the free sur- 
face. The results have shown that the calculations are fairly stable at 
relatively high Froude numbers. However, at relatively small Froude 
numbers, the Kelvin-Helmholtz instability sets in quickly and the free 
surface exhibits highly irregular shapes. It appears that either the 
integration procedures have to be refined or more sophisticated vor- 
tex sheets have to be used to calculate the small deformation of the 
free surface at small Froude numbers. 

The results obtained with a Froude number of Fv = 1.125 are 
shown in Figures 2.1 through 2.27. The most striking feature of these 
results is that the free surface rises rapidly to a height of about 1.25 
above the mean water level and captures the vortex pair. It has not 
been possible to carry out the calculations to larger times. The fact 
should be kept in mind that the point vortices used in the numerical 
calculations may have very little resemblance to the real vortices by 
the time they rise to the mean water level. In fact, the experiments 
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show that the vortices become quickly turbulent. There is at present 
no possibility of incorporating into a numerical analysis the motion of a 
turbulent vortex. 

Experiments were carried out at a maximum Froude number of 
0.35. The evolution of the free surface is shown in Figures 4.1 through 
4.12. The shape of the free surface at the time of maximum rise is 
plotted in Figure 5. Figures 4 and 5 show that the free surface rises to 
a maximum height of about 0.25. develops two strong scars, and then 
subsides as the vortices begin to move parallel to the free surface. The 
scar front is slightly ahead of the vortex center and reduces to a small 
depression at larger times. The results presented in Figures 4 and 5 
should form the basis of comparison for future numerical efforts at the 
corresponding Froude numbers. Additional calculations are underway 
with more sophisticated vortex sheets. The results presented herein 
are extremely encouraging and are expected to lead to a better under- 
standing of this extremely complex and challenging phenomenon. 
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VII. CONCLUSIONS 



The investigation reported herein warranted the following 
conclusions: 

1. The basic equations governing the motion of vortices in homo- 
geneous and density-stratified media have been developed and 
expressed in terms of dynamic and buoyant scaling laws; 

2. A novel experimental technique has been devised to generate a 
nearly two-dimensional vortex pair rising toward a free surface; 

3. Migration of a vortex pair toward a free surface gives rise to a 
bulge and two scars; 

4. The maximum rise at the free surface occurs when the vortex 
pair is at a depth of about 0.5; 

5. In the range of Froude numbers encountered in the experi- 
ments, the rise of the free surface is always followed by a fall as 
the vortices begin to move parallel to the free surface at large 
times; 

6. The numerical analysis of the free surface deformation through 
the use of discrete vortices leads to instabilities at low Froude 
numbers. These instabilities could have been removed with a 
suitable numerical filter or smoothing technique. However, the 
use of a relatively subjective smoothing procedure has been ruled 
out. The calculations at a Froude number of 1.125 have yielded 
fairly stable solutions and provided the basis for comparison with 
future experiments. 
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APPENDIX: FIGURES 
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and Boundary Conditions 
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Figure 2.1 Velocity Field for T* = -2.40 
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Figure 2.2 Velocity Field for T* = -2.30 
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Figure 2.3 Velocity Field for T* = -2.20 
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Figure 2.4 Velocity Field for T* = -2.10 
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Figure 2.5 Velocity Field for T* = -2.00 
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Figure 2.6 Velocity Field for T* = -1.90 
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Figure 2.7 Velocity Field for T* = -1.80 
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Figure 2.8 Velocity Field for T* = -1.70 
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Figure 2.9 Velocity Field for T* = -1.60 
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Figure 2.10 Velocity Field for T* = -1.50 
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Figure 2.1 1 Velocity Field for T* = -1.40 
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Figure 2.12 Velocity Field for T* = -1.30 
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Figure 2.13 Velocity Field for T* = -1.20 
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Figure 2.14 Velocity Field for T* = -1.10 
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Figure 2.15 Velocity Field for T* - -1.00 
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Figure 2.16 Velocity Field for T* = -0.90 
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Figure 2.17 Velocity Field for T* = -0.80 
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Figure 2.18 Velocity Field for T* = -0.70 
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Figure 2.19 Velocity Field for T* = -0.60 
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Figure 2.20 Velocity Field for T* = -0.50 
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Figure 2.21 Velocity Field for T* = -0.40 
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Figure 2.22 Velocity Field for T* = -0.30 
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Figure 2.23 Velocity Field for T* = -0.20 
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Figure 2.24 Velocity Field for T* = -0.10 
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Figure 2.25 Velocity Field for T* = 0.00 
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Figure 2.26 Velocity Field for T* = 0. 10 
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Figure 2.27 Velocity Field for T* = 0.20 



93 





li- 


\ 7 
/ \ 


! 

.i 

r- : 

i 

.1 

1 ! 

i 

l| 

3 


II* 


i\ /i 

II: | S r -+' 

i 1 % 

/ \ 


.U 

•i 

? **• 
* 

j 

;r 

1 

1 

i 

•>« 

i 

ij 

,1 


l;! 


I \ “7 

i x 
r 

|/ \ 


i 

«i 



I 

I 

t 

t 

t 



I 



> 

Ml 

> 

V- 

2 

* 

w 



r 



94 



Figure 3. Experimental Apparatus 
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Figure 4 . 1 Vortex Motion at T* = -1 .05 
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Figure 4.2 Vortex Motion at T* = -3.46 
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Figure 4.3 Vortex Motion at T* 
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Figure 4.4 Vortex Motion at T* = -2.73 
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Figure 4.5 Vortex Motion at I'* = -1.63 
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Figure 4.6 Vortex Motion at T* = -0.82 
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Figure 4.7 Vortex Motion atT* = -0.28 
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Figure 4.8 Vortex Motion atT* = -0.09 
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Figure 4.9 Vortex Motion at T* = 0.09 
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Figure 4.10 Vortex Motion at T* = 0.28 
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Figure 4.1 1 Vortex Motion at T* = 0.44 
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Figure 4.12 Vortex Motion at T* = 0.82 
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Figure 5, Free Surface Shape at Maximum Rise 
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